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SUMMARY 


A computational procedure has been developed which can be used 
to determine the flow properties in hypersonic helium wind tunnels 
in which real gas behavior is significant. In this procedure, a 
three-coefficient virial equation of state and the assumption of 
isentropic nozzle flow are employed to determine the tunnel 
reservoir, nozzle, throat, freestream and post-normal shock 
conditions. This method has been applied to a range of conditions 
which encompasses the operational capabilities of the NASA Langley 
22-Inch Mach 20 Helium Tunnel. Results are presented graphically 
in the form of real gas correction factors which can be applied to 
perfect gas calculations. Important thermodynamic properties of 
helium are also plotted versus pressure and temperature. The 
computational scheme used to determine the real -helium flow 
parameters has been incorporated into a FORTRAN code which is 
discussed in this paper. 

INTRODUCTION 

At high densities, intermolecular forces can have a 
significant effect on the pressure-density-temperature relationship 
and other thermodynamic properties of interest. For example, at 
typical reservoir conditions for the Langley 22— Inch Mach 20 Helium 
Tunnel of 137 atm. and 295° K, perfect gas theory yields a density 
of 22.65 kg/m 3 and a specific heat ratio of 1.667. However, the 
actual values of these properties at the specified conditions would 
be 21.26 kg/m 3 and 1.652, which amount to errors of +6.6% and +0.9% 
in the values obtained from perfect gas theory. For a measured 
test-section pitot pressure of 0.326 atm, the perfect gas Rayleigh 
pitot formula can be used to calculate a Mach number of 21.2, 
however the actual Mach number would be 21.6; thus, the perfect gas 
value is off by —1.8%. For a perfect gas the test section 
stagnation temperature would still be 295° K; however, the real gas 
temperature under these conditions would be 302° K. Thus, the 
perfect gas computation results in an error of +2.3%. This tunnel 
can be operated at much more extreme conditions than these, which 
would cause a much greater deviation from perfect gas behavior. 

A method has been devised which accounts for real gas effects 
through the use of the virial form of the equation of state and 
real gas thermodynamic relationships for expanding flows. In this 
method, the standard assumption of isentropic flow is made. It is 
also assumed that the flow is remains free of condensation during 
the expansion process, that it remains at temperatures below that 
at which electronic excitation becomes a significant factor 
(T < 10000° K) 1 , and that the helium density remains below the 

critical point (p CR = 69.64 kg/m 3 ) 2 . These conditions are all 

satisfied in the normal operating range of Langley's Mach 20 helium 
wind tunnel. 

Analysis of real— helium flow parameters has been carried out 
by Erickson 3,4 at moderate enthalpies and at high enthalpies by 
Miller and Wilder 5 . In Erickson's work, the Beattie-Bridgeman 
equation of state was used to obtain results up to 410 atm. and 
590° K. This work was later extended to obtain approximate results 
up to 3400 atm. and 5800° K. Miller and Wilder were concerned with 


1 


the prediction of hypervelocity flow parameters for Hotshot tunnels 
and used a virial equation of state with coefficients determined by 
theoretical analysis of high-temperature gas behavior to obtain 
results for up to 3600 atm. and 15,000’K. In the present study, a 
virial equation of state is also used, but with virial coefficients 
derived from both analysis of low and moderate temperature 
experimental data and high-temperature theoretical results. This 

? r °v ldeS accurate results both in the operational range of 
the 22- In ch Mach 20 Helium Tunnel (20 - 285 atm, 270 - 600° K) 6 and 
in the high pressure, high temperature range of computations 
- rrie in reference (5). However, in this paper the emphasis 
gJeate^thln^O ° f 1 ~ 400 atTa * and 50 “ 6 °0° K with Mach numbers 


SYMBOLS 


a 

B 

C 


D 


e 

h 

M 

P 

R 

Re 

s 

U 

W He 

Z 


Y 

P 


speed of sound (m/s) 

second virial coefficient (m 3 /kg) 

third virial coefficient (m 3 /kg) 2 

specific heat at constant pressure (J/kg-'K) 

specific heat at constant volume (J/kg-°K) 

fourth virial coefficient (m 3 /kg) 3 

internal energy (J/kg) 

enthalpy (J/kg) 

Mach number 
pressure (N/m 2 ) 

gas constant for helium (2077.1 J/kg-°K) 
Reynolds number 
entropy (J/kg) 
velocity (m/s) 

molecular weight of helium (4.0026 g/g-mole) 

compressibility factor 

specific heat ratio 

coefficient of viscosity (kg/m-s) 

density (kg/m 3 ) 
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Subscripts : 

0 total 

1 freestream conditions 

2 post-shock static conditions 

P perfect gas value 

r reference value 

t,l reservoir stagnation conditions 

t,2 post-shock stagnation conditions 

x conditions at arbitrary point 


Superscripts : 

* throat conditions 

VIRIAL EQUATION OF STATE 

Flow properties are determined through the use of the virial 
equation of state, real gas relationships for entropy, enthalpy, 
specific heats and speed of sound, and the shock conservation 
relations for mass, momentum and energy. 

The equation used to represent the pressure-density- 
temperature relationship of real helium gas is the virial equation 
of state 


p = ZpRT = pRT[l + p B(T) + p 2 C(T) + p 3 D(T) + . . .] (i) 


A relatively large amount of experimental data 7 ' 11 on the second 
virial coefficient, B(T) , exists for low and moderate temperatures 
(up to about 1000° K) . Experimental results for the third virial 
coeff icient 7,10,11 , C(T) , are few and show a great deal of scatter. 
No experimental data on the fourth and higher virial coefficients 
are available. At high temperatures (1000° K+) , theoretical 
results for the second, third and fourth virial coefficients were 
derived from models of the intermolecular force potential of helium 
presented by Amdur and Mason 12 . 

Expressions were derived for B(T) and C(T) which match both 
the theoretical results and the experimental data, while that for 
D(T) was based is a theoretical relationship taken from reference 
(12) . At the conditions dealt with in this work, the contribution 
of C(T) to the equation of state is relatively small and that of 
D(T) essentially negligible, thus the virial equation of state is 
truncated at the third virial coefficient. If calculations are to 
be carried out at high temperatures and densities, the fourth 
virial coefficient should be included. 
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The expressions derived for the virial coefficients are 

B(T) = b 0 + + b 2 (T~ 3/i ) + b 3 (T- 5/i ) + b 4 (T- 7 ' 4 ) (2) 

T Z 1300°Jir 
b 0 = -13.4067 

b 3 = 165.4459 

b 2 = -1357.92 

b 3 = 5959.061 
b 4 = -12340.8 


r>i3oo°ic 
b Q = 1.178236 
b 3 = -7 .57134 
b 2 = 5225.701 
b 3 = -188923. 
i> 4 = 2460461. 


c{t) = c 0 + c^r- 174 ) + c 2 (T- 3/i ) + c 3 (r- 5/4 ) 

c 0 = -13.7898 
c x = 139.7339 
c 2 = 8114.259 
c 3 - -17456.9 


where the units of T, B(T) and C(T) are °K, cm 3 /mole and (cni 3 /mole) 2 
respectively. Expressions for B(T) , C(T) and D(T) given by 
previous authors are: 

Amdur and Mason 12 


B(T ) = 1.3436 x 10- 2 (15.8922 - InT) 3 (4) 

C(T) = 9.0263 x 10~ 5 ( 15 . 8922 - InT) 6 (5) 

D(T) = 7.0341 x 10" 7 ( 15 . 8922 - InT) 9 (6) 


4 



(V) 


Harrison 13 

B(T) = 1.3436 x 10 * 2 (15 . 8922 - InT) 3 - 4 . 39exp [-2 . 4177 x 10' 3 T] 


Miller and Wilder 5 

B{T) = 1.3436 x 10 ~ 2 (15. 8922 - InT) 3 - 8 . 04exp [ -3 . 7 156 x 10~ 3 T] 


( 8 ) 

The values of B (T) and C(T) computed from equations (2) and (3) are 
plotted in figures (1) and (2). Also plotted are the experimental 
data and theoretical results on which these expressions are based 
and the values of B (T) and C(T) obtained using equations (4), (5), 

(7) and (8) . Note that for ease of comparision with existing data, 
the virial coefficients are plotted in terms of (cm 3 /mole) , 
where: (cm 3 /mole) = W He /1000 (m 3 /kg) . 

At high temperatures, the various expressions for B(T) 
approach the same values because they are all based on Amdur and 
Mason's theoretical computations. However, in that work only the 
effect of the intermolecular repulsive force, which dominates at 
high temperature, is taken into account; the effect of the 
intermolecular attractive force, which dominates at low 
temperature, is ignored. Thus, the values of B(T) from equation 
(5) become increasingly inaccurate at lower temperatures. Amdur 
and Mason's expression for B(T) was later corrected empirically to 
reflect the effect of the attractive force through the addition of 
logarithmic terms as shown in equations (7) and (8) . 

Between 500 and 1000° K, the expression used herein for B (T) 
and those from the other references can be used to obtain 
approximately the same values. However, as shown in figure (1), 
these other expressions cannot be used to model the behavior of 
B(T) below 200°. Between 1000 and 2000° K, the results obtained 
from these expressions also differ significantly from the 
experimental data. The effects of the second (and higher) virial 
coefficients are small in the first range of temperatures stated 
for the densities generally associated with them in hypersonic 
flows, and thus their accuracy is not of great concern. However, 
for temperatures between 1000 and 2000° K the influence of B(T) is 
not negligible, and the use of these expressions may introduce a 
considerable error. The expression presented in this work, 
equation (2) , can be used to accurately calculate the value of B(T) 
across the temperature range of interest. 

As seen in figure (2) , third virial coefficient data shows a 
considerable amount of scatter, but the accuracy of the expression 
for C(T) is not as great of a concern since its effects on the 
virial equation of state are small over the temperature range of 
interest. The expression used in the present study, equation (3), 
is fitted to the high temperature data of Amdur and Mason and also 
agrees with the trend of the experimental data at low and moderate 
temperatures . 

The accuracy of these coefficients was checked by comparing 
the values of the resulting compressibility coefficient, Z, with 
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those given in the tabulation of helium properties prepared by the 
International Union of Pure and Applied Chemistry* (IUPAC) . In 
figures (3a) and (3b) the data are compared along nine isobars 
ranging from 20 MPa to 0.001 MPa (197 atm to 0.00987 atm). The 
compressibility coefficient is also plotted versus pressure and 
temperature in figure (4) . The only region where the computed 
values deviate significantly from the IUPAC values is in the very 
low temperature-very high density range, and this is generally not 
a condition dealt with in hypersonic research. 

THERMODYNAMIC PROPERTIES 

Real gas thermodynamic relationships based on the virial form 
of the equation of state were derived in references (1) and (3). 
The derivations are repeated here for completeness. 

Enthalpy 

Expressing enthalpy as a function of density and temperature gives 



And by definition 

h = e + (10) 

P 


Differentiating this equation with respect to T at constant p 



( 11 ) 


And by definition 





( 12 ) 


By the first law of thermodynamics 

6 q + 6 w = de 


( 13 ) 


Or for isentropic compression or expansion 
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( 14 ) 


Tds + — pdp = de 

P 2 


Differentiating equation (10) and substituting into equation (14) 
yields 


dh = Tds + —dp 

P 


(15) 


At constant temperature, equations (9) and (15) give 

Up l T \ 5p/ r P\ 6 P/r 


(16) 


And by the Maxwell relation 


I ds \ _ -1/ 6p\ 

Up I t p 2 U t} p 


(17) 


Then 


/ _ 

j^IAp\ 

+ i/8p.\ 

Up / r 

P 2 UrJ p 

p Up / 


(18) 


By substitution of equations (11) and (18), equation (9) becomes 


dh = 


^ .±I*E\ ] 

dT + 

— 1 

i£) * — 

?>P\ ' 

Q.’ 

«o 

Q_ 

> 

i 


p 2 ' 

i8r|, P 1 

6 P / T 


dp 


(19) 


Where, from equation (1) , the virial equation of state 


(H). ■ 


pi? 


1 + p(b(D + r dB d y } | + p 2 {c(T) + r ) 


( 20 ) 


and 


= i?T[l+2pS(r) + 3 p 2 C( T) ] (21) 

°P I T 
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Since enthalpy is a potential function, the integral of dh is 
independent of the path of integration. Therefore, integrate 
equation (19) in two steps: the first step at constant density and 
the second at constant temperature. 


f T ' P dh = f T ' Pl dh + f T ' P dh 

J T I .p r JT. p r 


( 22 ) 


For the first step, integrate to the desired temperature at an 
arbitrary, constant density which is low enough that perfect gas 
behavior can be assumed. Thus, for this step, equation (11) 
reduces to: 



c v + R 


= c r 


( 11 ) 


then 



c p (T - T x ) 



(23) 


The temperature term evaluated at the reference temperature is 
eliminated by the value of enthalpy at the reference condition 
(h = 1548.2 kJ/kg at T = 298.15 °K) 2 since it is chosen to be a 
point of perfect gas behavior. 

For the second step, integrate from the reference density to 
the desired density at constant temperature to get 



L 


T. P 
T.Pr 


Tl(*P 

bT‘ 


-4<4£>P - f 

2r dC(T) 


+ RT-P- 

2 


2 C{T) - T 


dp = RTp 

P 


B(T) 


dT 


Pr 


T dB[T) 

dT 


p 

-Pr 


(24) 


At the reference density, perfect gas behavior is assumed. Thus 
the virial coefficients and their derivatives at this point can be 
neglected. Then, by substitution of equations (23) and (24) into 
(22) , the expression for enthalpy becomes 


h 



p[B(T) 


rp dB ( T) \ 

dT I 


-^-l2C(T) - T \ 

2 \ dT ) 


(25) 


Entropy 

An equation for entropy can also be derived by first 
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(26) 


expressing s as a function of p and T 


ds 




dp 


The second term has already been given by equation (17) . 
internal energy as a function of p and T gives 


de 




dp 


Defining 


(27) 


Substituting equations (12) and (14) gives 


ds 


-^dr + 

T 



(28) 


And since equations (26) and (28) are equivalent 

6s \ _ 

6 Tj p T 


(29) 


By substituting equations (17) and (29) into equation (26) , the 
expression for entropy is obtained 


ds 


-^dT - 

T 




(30) 


The expression for entropy can be integrated in the same manner as 
that for enthalpy to get 


s = R\ 


In ( T) 


- lnp - p| 


B ( T) + T 


dB ( T) 

dr 


i( c{T) 


+ T- 


dC(T) 

dT 


+ s 


xef 


(31) 


The value of the constant, s ref , is determined by substituting 
the values for reference entropy, temperature and density (31.489 
kJ/kg-°K at 298.15° K, 0.16361 kg/m 3 , respectively) 2 into equation 
(31) to get s ref = 9.977 kJ/kg-°K. 

Other Properties 

The specific heat at constant volume can then be found using 
equations (29) and (31) 
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cv = 


Yp - 1 


- RT 


pf 2 -^in + T d 2 B{T) ) P 2 f dC{T) + T d 2 C{T) ) 

( dT 2 ) 2 f dr dr 2 J. 


(32) 

The speed of sound and specific heat at constant pressure are found 
through the thermodynamic relationships 


= IAp\ = c p/ &P\ 

Up/ a c„\5p] 


(33) 


Up 


~T ibp\ 


c v p 2 \ 6P/ p 

UW P 

\ = JAp\ 

/H\ 

It U Tf 

pUp/p 


(34) 


(35) 


Substitution of equations (34) and (35) into (33) leads to 


(- f£\ + 

t /Spy 

Up It 

c/Ur|, 


(36) 


The speed of sound can then be computed using equations (20) and 
(21) with equation (3 6) . 

The specific heat at constant pressure is found by 
substituting equation (35) into (34) to obtain 



(37) 


The final relations required are the one-dimensional shock 
conservation equations 


mass: 


P2^2 Pl^l 


(38) 


momentum : 


Pi + Pl^l = P 2 + P 2 U 2 2 


(39) 
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energy : 


(40) 


ul u 2 2 

= h x + — ^h 2 + — 


In figure (5) , values of enthalpy and entropy calculated using 
the procedure described are plotted along lines of constant 
pressure and temperature. The speed of sound is plotted in figure 
(6) and the specific heats in figures (7) and (8) . 

COMPUTATIONAL PROCEDURE 

The relationships from the previous section can be used to 
determine the reservoir, nozzle, throat, freestream and post-normal 
shock flow parameters in a helium wind tunnel from the measured 
reservoir pressure and temperature and the measured pitot pressure 
or freestream Mach number. The procedure described herein has been 
incorporated into a FORTRAN program called HEPROPS which is listed 
in the appendix. 

Reservoir Conditions 

Given the reservoir pressure and temperature, which are the 
properties normally measured during tunnel operation, and assuming 
zero flow velocity in the reservoir, the corresponding reservoir 
density can be obtained through Newtonian iteration on the virial 
equation of state, where 

F = p 0 n -p 0 i?r 0 [i + p 0 b(t 0 ) + poC(r 0 )] 


F'=-RT 0 [l + 2p 0 B(T 0 ) + 3poC(T 0 )] 


At this point, if the freestream Mach number is known, the 
computations can be carried out in a straightforward manner for the 
throat, freestream, and post-normal shock conditions. If instead, 
the pitot pressure is given, a Mach number must be assumed and all 
computations carried out for that Mach number. Newtonian iteration 
must then be carried out on the Mach number until the measured and 
computed pitot pressures are in agreement. 

Conditions at An Arbitrary Point 

Conditions at any point in the nozzle can be determined by 
specifying a Mach number and assuming that the expansion through 
the nozzle is an isentropic. This procedure is used to determine 
the conditions at the throat (M* = 1) and in the freestream (M, 


(41) 

(42) 

(43) 
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either given or assumed during iteration) . Additionally, the 
option exists in the HEPROPS program to compute the conditions at 
any arbitrary point (denoted "x") in the nozzle by specifying a 
Mach number for that point. 

From the total enthalpy and entropy, h„ and s 0 , (which can be 
computed from the supply conditions) and a specified Mach number 
(for instance, at the throat, M x = 1 ) . Static density and 
temperature, p x and T, can then be determined through Newtonian 
iteration on two— variables, where 


G = 


Px 


( 44 ) 


W 


(h Q - RT X 


Yd - 1 


+ Px 


B(rj - 

x dr„ 


M 




2 C ( T ) - T x - °^ Tx) ' 

X x dT 

X / 


S Q - R 


TT - ln P* - P.ffilr,) * T^IA.y^c( T x )T x ^SlA^ 


W a * 1 = w 1 


fityi - 1 
6 <J 


( 45 ) 

( 46 ) 


Where the Jacobian of W is given by 





bw 2 

6 tr\ 



Sg 1 

8 G 


8 w 1 

8w 2 



bg 2 

bg 2 


( 47 ) 


Post— Shock Static Conditions 

Post— shock static conditions can be determined from the 
computed freestream conditions using the shock conservation 
eguations. From these equations, implicit expressions for p, and 
T2 can be derived and then solved iteratively, where c 



( 48 ) 
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w 


Oh - RT 2 


u| 


‘l-E— + p Jb(T 2 ) ~ T 2 dB{T A 


Yp - 1 


dr. 


i * - ! ( 


p| 2 C{T 2 ) - T 2 


dC(T 2 ) 

dT 0 


p. + p„ui - p 2 i?r 2 [i + p 2 B(T 2 ) + p2C(T 2 )] - 


(p.u°°) 2 


( 49 ) 


Post-Shock Stagnation Conditions 

By assuming isentropic compression from the post-shock static 
conditions, p t 2 and T t 2 can be computed using 


(ho ~ RT ; 


t,2 


Y / v dB{T t2 ) 

+Pt , 2 B(T t , 2 ) 




/ dC(T t 2 ) \ . 

- Pu 2 [2C(T t , 2 ) - r c , 2 dr ^ 2 -j ) 


(s, - R 


Yp - 1 


- ln(p t/2 ) - pt, 2 


b(t C i2 ) + r 


dfl(T c>2 ) 


dr, 


t,2 


P t, 2 


Z , % dC(T t 2 ) \ 


(50) 


a = 


P t, 2 

r t .2. 


(51) 


At each point, the remaining thermodynamic properties can be 
evaluated using the eguations from the previous section and by 

M = — (52) 

a 


T-, P ux 
Re = -B- 

P 


(53) 
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from reference (14) 


T <; i.2*y 

H = (2.1630 - 26 .665 T + 120.54 T 2 - 187.41 T 3 
+ 126 .82 T 4 - 31 . 823 T 5 ) xlO -7 - jcgr 

m-s 


1.2*K < T < 3 

JA = (5.02 - 3.2241 T + 2.0308 T 2 - 0.22351 T 2 ) xlO -7 

m-s 


3 ■ 6 "at < r < 8°y 

n = (-1.5691 + 3.4167 T - 0.10317 T 2 ) xlO' 7 -iS 

m-s 


T . 1 Bl K 

H = 5 . 02 3x1 O' 7 ( 54 ) 

m-s 


The expression used for viscosity 14 is plotted versus 
temperature in figures (9a) and (9b) and is compared with 
experimental and theoretical data 7 - 12 ' 1 *' 18 . 

At this point, the computation is complete if the Mach number 
was originally specified. However, if instead the pitot pressure 
was given, Newtonian iteration for the Mach number must be carried 
out where 


F = p 


t, 2 


calc 


Pt.i 


(55) 


f' is evaluated numerically 

f>P t 2 , 

L c calc 


F = 


Sal 


Pt,2^» + AMJ ~ P t , 2 (Mj 
Am 


(56) 


bC + 1 = K? - JL 

F' 


(57) 


It should be noted that this method provides a most accurate 
means of computing flow properties over a wide range of conditions 
since the assumption made in previous works of perfect gas behavior 
in the freestream (which is fairly accurate for the Langley 22-Inch 
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Mach 20 Helium Tunnel) is not made in these computations. 

RESULTS AMD DISCUSSION 

The HEPROPS program was used to determine helium flow 
properties for reservoir conditions of 1 to 400 atm. pressure and 
50 to 600 °K temperature, which encompasses the operating range of 
the Langley 22-Inch Mach 20 Helium Tunnel (20 - 285 atm. and 290 - 
600' K) . Results are tabulated graphically in figures (10) through 
(13) . In each figure are shown correction factors for ratios of 
thermodynamic properties at different points in the flow. These 
correction factors represent the actual ratio for the given 
thermodynamic property (as determined using HEPROPS) divided by the 
ratio calculated from perfect gas theory. The correction factors 
are given for Mach 10 and 20. The correction factors for Mach 
numbers greater than 20 are essentially constant and can be taken 
as their values at Mach 20. These figures should not be used to 
estimate correction factors for flows with Mach numbers much less 
than 10; however the operating conditions of low Mach number flows 
are generally such that real gas deviations are insignificant. 

As part of the computational procedure, the specific heats at 
constant pressure and temperature, the speed of sound, and the 
entropy and enthalpy were calculated across the entire range of 
flow conditions. The values of the properties can be estimated at 
any given pressure and temperature using figures (5) to (8) . Also, 
the virial coefficients and the compressibility coefficient are 
given in figures (1) through (4) . Viscosity is given in figures 
(9a) and (9b) . 

Ideally, the HEPROPS program should be used to evaluate wind 
tunnel flow properties. However, if desired, flow properties can 
be estimated by applying these real gas correction factors to 
perfect gas calculations. As an example, consider a helium wind 
tunnel flow for which the freest ream Mach number, M.,, is equal to 
20, and the measured supply pressure and temperature are p 0 = 30.4 
MPa (300 atm) and T 0 = 300 K. The supply density is 

= P ° 

° Z 0 RT 0 


From figure (4) 


Po 


3 0,4 MPa 

(1.14) -(2079 J/kg-K) • (300 K) 


42 . 8 kg/m 2 


The freestream pressure is 






y P 


1 
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From figure (11a) 


Pi = 


(1.14) 


3 0.4 MPa 


1 + 


5/3 - 1 


5/3 


20 


21 5/3 


= 166 Pa 


Note from figure (lie) that the correction factor for the 
freestream Mach number was found to be between 0.99 and 1 through 
most of the range of conditions studied, thus for rough estimates 
no correction to the Mach number is required. 

The post-normal shock stagnation temperature is 



From figure (13b) 


T t 2 = (1.06) - 00K = 318 K 
1 


CONCLUDING REMARKS 

A method for computing real gas flow parameters in helium wind 
tunnels has been developed which utilizes a three-coefficient 
equation of state. Comparisons of the compressibility 
factors obtained using this equation were made with those in 
published tabulations and were found to be in very close agreement. 
This method has been incorporated into a FORTRAN program designed 
to be used in the determination of wind tunnel flow properties in 
the reservoir , nozzle, throat, freestream and post-normal shock 
This program can also be used to calculate thermodynamic 
properties at any arbitrary point in the nozzle which is specified 
by a Mach number. Correction factors have been tabulated for the 
range of conditions in which the Langley 22-Inch Mach 20 Helium 
Tunnel is operated. These correction factors can be applied to 
perfect gas theory computations to estimate real flow properties. 
Values for relevant thermodynamic properties of helium as functions 
of pressure and temperature are also presented. 
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Figure 2. — Third Virial Coefficient as a Function of Temperature 
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Temperature, T, (K) 

Figure 3. - continued 

(b) P * 2.5 MPa 
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Figure 4. - Compressibility Coefficient vs Temperature and Pressure 
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Figure 5. — Entropy vs Enthalpy Chart 
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Figure 6. - Speed of Sound as a Function of Temperature and Pressure 
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Figure 7. - Specific Heat at Constant Pressure 







Figure 8. - Specific Heot at Constant Volume 
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Figure 9. - Viscosity as a Function of Temperature 
(a) High Temperature Range 
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Figure 10. - Ratios of Real to Ideal Flow Parameters 
(a) Pressure Ratio, (po/p*) 
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Supply Pressure, p 0f (atm) 

2 10. — Continued 

(b) Temperature Ratio. (T 0 /T*) 
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Supply Pressure, p 0 , (atm) 

Figure 10. - Continued 

(c) Density Ratio, (rho 0 /rho*) 
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Supply Pressure, p 0 , (atm) 

Figure 10. - Continued 

(d) Speed of Sound Ratio (a 0 /a*) 
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100 150 200 250 300 350 

Pressure, p 0 , (atm) 

Figure 10. — Continued 

(e) Mach Number Ratio, (M* p /M*) 
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100 150 200 250 300 350 400 

Supply Pressure, p 0 , (atm) 

Figure 11. — Continued 

(b) Temperature Ratio, (T 0 /Ti) 
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Supply Pressure, p 0 , (atm) 

Figure 11. - Continued 

(d) Speed of Sound Ratio (ao/a,) 
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Figure 12. - Ratios of Real to Ideal Flow Parameters 
(a) Pressure Ratio, (P 0 /P 2 ) 
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Supply Pressure, p D> (atm) 

Figure 12. - Continued 

(b) Temperature Ratio, (T 0 /T 2 ) 
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Supply Pressure, p D , (atm) 

Figure 12. — Continued 

(d) Speed of Sound Ratio, (a 0 /a 2 ) 
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Figure 13. - Ratios of Real to Ideal Flow Parameters 
(a) Pressure Ratio, (pd/pu) 
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Figure 13. - Continued 

(b) Temperature Ratio, (To/T u ) 
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100 150 200 250 300 350 

Supply Pressure, p 0 , (atm) 

Figure 13. — Continued 

(c) Density Ratio, (rho 0 /rhot, 2 ) 


APPENDIX 


PROGRAM HEPROPS 
C BRIAN R. HOLLIS 

C MARCH 1992 

C — — — 

C - PROGRAM HEPROPS IS DESIGNED TO EVALUATE THE WIND 

C = TUNNEL FLOW PROPERTIES OF HELIUM GAS. IT IS 

C = ASSUMED THAT THE EXPANSION OF THE TEST GAS THROUGH THE 

C = TUNNEL IS ISENTROPIC AND THAT TOTAL ENTHALPY IS CONSERVED. = 

C - THE HELIUM IS TREATED AS A REAL GAS THROUGHOUT THE 

C - EXPANSION. THIS PROGRAM CAN ALSO BE USED TO EVALUATE THE 

C - PROPERTIES OF REAL HELIUM AT A GIVEN TEMPERATURE AND 

C - PRESSURE OR TO COMPUTE FLOW PROPERTIES AT AN 

C - ARBITRARY MACH NUMBER IN THE TUNNEL EXPANSION SECTION 

C - 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


“ - PROGRAM VARIABLES AND CONSTANTS - 

A - SPEED OF SOUND 

CP - SPECIFIC HEAT AT CONSTANT PRESSURE 

CV - SPECIFIC HEAT AT CONSTANT VOLUME 

- flag = SETS OPERATIONS TO BE CARRIED OUT IN PROGRAM 
= HREF - REFERENCE ENTHALPY AT 298.15 K, latm 

K - SPECIFIC HEAT RATIO 

- MACH - MACH NUMBER 

MHE - MOLAR MASS OF HELIUM 
MU - VISCOSITY 

P - PRESSURE 

- PITOT - TEST SECTION PITOT PRESSURE 

R - UNIVERSAL GAS CONSTANT 
RHO - DENSITY 

- SREF - REFERENCE ENTROPY AT 298.15 K, 1 atm 

VB - SECOND VIRIAL COEFFICIENT 

“ VBP ” FIRST DERIVATIVE OF VB WITH RESPECT TO TEMPERATURE = 
= VBPP - SECOND DERIVATIVE OF VB WITH RESPECT TO TEMPERATURE = 
” VB3P - THIRD DERIVATIVE OF VB WITH RESPECT TO TEMPERATURE = 

- VC - THIRD VIRIAL COEFFICIENT 

- VCP - FIRST DERIVATIVE OF VC WITH RESPECT TO TEMPERATURE = 
= VCPP - SECOND DERIVATIVE OF VC WITH RESPECT TO TEMPERATURE = 
= VC3P - THIRD DERIVATIVE OF VC WITH RESPECT TO TEMPERATURE = 

Z - COMPRESSIBILTY FACTOR 

= THE VARIOUS ARRAYS APPEARING IN THE PROGRAM ARE FOR = 

STORAGE OF DATA DURING ITERATIONS. 

- VARIABLE SUFFIXES - 

0 - SUPPLY RESERVOIR 

1 = FREESTREAM 

2 - POST NORMAL SHOCK 
S - SONIC POINT 

X - ARBITRARY POINT IN NOZZLE SPECIFIED BY MACH NUMBER = 
T - STAGNATION CONDITION 
P - EQUIVALENT PERFECT GAS CONDITION 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


THUS, P2 IS STATIC PRESSURE BEHIND THE NORMAL SHOCK = 
WHILE PT2P IS THE EQUIVALENT PERFECT GAS PITOT PRESSURE . = 
CP, CV, AND K WITHOUT SUBSCRIPTS ARE THE PERFECT GAS 
VALUES FOR HELIUM. 


- UNITS - 

ALL COMPUTATIONS IN THE PROGRAM ARE CARRIED OUT IN = 
S.I. UNITS. INPUT DATA UNITS ARE NOTED BELOW. THE 
UNITS OF THE OUTPUT ARE UP TO THE USER TO FORMAT. 


IMPLICIT REAL (A-Z) 

INTEGER I , J , L,N 

DIMENSION W(2) ,G(2) ,DWDG(4,4) , OLDW(2) , INV(2 , 2) , OLDG ( 2 ) , FP ( 2 ) 
COMMON K ,MHE ,R , SREF 

***** CONSTANTS ***** 

HREF = 0 
MHE = 4.0026 
R = 8314 
SREF « 9.977E3 


***** IDEAL GAS VALUES ***** 
K = 5/3. 

CP = K*(R/MHE)/(K-1) 

CV - (R/MHE)/ (K- 1) 


************************************************* 
* - DATA INPUT - * 


* FLAG = 0 COMPUTE PROPERTIES USING MACHl * 

* =1 COMPUTE PROPERTIES USING PITOT * 

* =2 JUST COMPUTE FLOW PROPERTIES AT * 

* AN ARBITRARY, MACH NUMBER SPECIFIED * 

* POINT IN THE NOZZLE. * 

* =3 OUTPUT ALL CORRECTION FACTORS * 

* FOR A GIVEN MACH NUMBER. * 

* MACHl - TEST SECTION MACH NUMBER * 

* P01 = SUPPLY PRESSURE (atm) * 

* PITOT = TEST SECTION PITOT PRESSURE (atm) * 

* T01 - SUPPLY TEMPERATURE (R) * 


* ALL QUANTITIES ARE CONVERTED TO S.I. UNITS * 
************************************************* 


*********************************************** 

* USERS CAN SUPPLY A DATA INPUT ROUTINE * 

* APPLICABLE TO THEIR SYSTEM AT THIS POINT, * 

* THE ROUTINE MUST READ IN THE QUANTITIES * 

* NOTED ABOVE AND CONVERT TO THEM S.I. UNITS. * 
*********************************************** 


***** SAMPLE INPUT FOR TESTING PROGRAM ***** 
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PO - 300*1. 01325E5 
TO - 300 
FLAG - 3 
MACH1 = 20 

PITOT = 1.01325E5*6. 04042/14. 697 


C 

C 

c 

c 

c 

c 

c 

c 


- SUPPLY CONDITIONS - 

SUPPLY RESERVOIR HELIUM BEHAVES AS A REAL GAS 
SUPPLY RESERVOIR PRESSURE AND TEMPERATURE 
ARE MEASURED. 


C ************************************************ 

C * SOLVE FOR DENSITY FROM VIRIAL GAS EQUATION * 

C * OF STATE. USE IDEAL GAS EQUATION OF STATE * 

C * FOR INITIAL ESTIMATE. * 

C ************************************************ 

RHOO - P0/( ( R/MHE )*T0) 

CALL VIRIAL (TO , VBO , VBPO , VBPPO , VB3P0 , VCO , VCPO , VCPPO , VC3P0) 

100 F - PO - RHOO* ( R/MHE ) *T0* ( 1 + VBO*RHOO + VCO*RHOO**2) 

DFDRHO - - (R/MHE) *T0*(1 + 2*VBO*RHOO + 3*VCO*RHOO**2) 

OLDRHO - RHOO 

RHOO = RHOO - F/DFDRHO 

ERR = ABS(( RHOO -OLDRHO) /OLDRHO) 

IF (ERR.GT.O. IE-5) GOTO 100 

C ***** COMPUTE RESERVOIR THERMODYNAMIC PROPERTIES ***** 

ZO = 1 + VBO*RHOO + VCO*RHOO**2 

HO = (R/MHE) *T0*(K/ (K-l) + RHOO*(VBO - T0*VBP0) + RHOO**2 
C /2*(2*VC0 - T0*VCP0) ) + HREF 

SO - (R/MHE) *(1/ (K-1)*LOG(TO) - LOG(RHOO) - RHOO*(VBO + TO* 

C VBPO) - RHOO**2/2*(VCO + TO*VCPO)) + SREF 

CVO = (R/MHE)/ (K-l) - ( R/MHE ) *T0* (RHOO* ( 2*VBP0 + T0*VBPP0) + 

C RHOO**2/2*(2*VCPO + T0*VCPP0)) 

DPDT = RHOO*(R/MHE)*(1 + RHOO* (VBO + TO*VBPO) + RH00**2*(VC0 + 
C TO*VCPO) ) 

DPDRHO - ( R/MHE ) *T0* ( 1 + 2*RHOO*VBO + 3*RHOO**2*VCO) 

CPO = CVO + (TO/RHOO**2)*DPDT**2/DPDRHO 
KO - CPO/CVO 

AO - (CPO* (DPDRHO) /CVO) **.5 
IF ( FLAG . EQ . 2 ) GOTO 750 

C *********************************************** 

C * USERS 'S CHOICE OF INPUTTING FREESTREAM MACH * 

C * NUMBER OR PITOT PRESSURE. GIVEN PITOT * 

C * PRESSURE, MUST ITERATE FOR MACH NUMBER, * 

C * GIVEN MACH NUMBER, NO ITERATION REQUIRED. * 

C *********************************************** 
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n o 


500 CONTINUE 

L = 1 

DELMACH = MACH1/1E5 
IF ( FLAG . EQ . 1 ) THEN 

MACH1 = MACH1 + DELMACH 
END IF 

510 CONTINUE 


C 

c 

c 

c 

c 

c 


- FREESTREAM CONDITIONS - 

ASSUME FREESTREAM BEHAVIOR IS REAL, TOTAL 
ENTHALPY AND ENTROPY ARE CONSERVED. 


SI = SO 
HTl = HO 

C ***** COMPUTE REAL GAS PROPERTIES ***** 

C ***** INITIAL GUESS AT (RH01, Tl) FROM PERFECT GAS EQUATIONS ***** 

Tl = T0/(1 + (K-1)/2*MACH1**2) 

RH01 = RHOO/(1 + (K- 1) /2*MACH1**2)**( 1/ (K- 1) ) 

CALL REALEXP ( Pi , Tl , RH01 , Z1 , MACH1 , U1 , Al , HI , HTl ,Sl, RE1 , MUl , Q1 ) 


C 

C 

c 

c 

c 

c 

c 


- POST -SHOCK STATIC CONDITIONS - 

ASSUME POST- SHOCK STATIC GAS BEHAVIOR 
IS REAL AND TOTAL ENTHALPY IS CONSERVED. 


HT2 = HO 


C ******************************************************* 

C * ITERATE FOR RH02 , T2 TO SATISFY CONSERVATION OF * 

C * MOMENTUM ACROSS SHOCK AND TOTAL ENTHALPY USING TWO- * 

C * VARIABLE (RH02 , T2) SHOOTING TECHNIQUE. INTIAL * 

C * ESTIMATE FROM PERFECT GAS RELATIONS. * 

C ******************************************************* 

RH02 = RH01*(K+l)*MACHl**2/( (K-1)*MACH1**2 + 2) 

T2 = Tl* ( 2*K*MACH1**2 - (K-l) )*( (K-1)*MACH1**2 + 2)/((K+l) 
C **2*MACH1**2 ) 

G(l) = RH02 
G(2) = T2 
150 CONTINUE 

CALL VIRIAL (T2 , VB2 ,VBP2 , VBPP2 ,VB3P2 , VC2 ,VCP2 ,VCPP2 ,VC3P2) 
W(l) - PI + RH01*U1**2 - RH02* ( R/MHE ) *T2* ( 1 + RH02*VB2 + 

C rh 02**2*VC2) - (RH0l*Ul)**2/RH02 
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W(2) - HT2 - (R/MHE)*T2*(K/ (K-l) + RH02*(VB2 - T2*VBP2) + 

C RH02**2/2*(2*VC2 - T2*VCP2)) - (RH01*Ul/RH02)**2/2 

C ***** COMPUTE JACOBIAN OF MATRIX [W] ***** 

DWDG (1,1) = - (R/MHE)*T2*(1 + 2*RH02*VB2 + 3*RH02**2*VC2 ) + 

C (RH01*U1/RH02 ) **2 

DWDG(1,2) = - (R/MHE)*(RH02 + RH02**2*(VB2 + T2*VBP2) + RH02**3 
C *(VC2 + T2*VCP2) ) 

DWDG (2,1) = - (R/MHE)*T2*((VB2 - T2*VBP2) + RH02*(2*VC2 - T2* 

C VCP2 ) ) + (RH01*U1)**2/RH02**3 

DWDG (2, 2) = - (R/MHE)*(K/(K-1) + RH02*(VB2 - T2*VBP2 + T2**2 

C *VBPP2) + RH02**2/2*(VC2 - T2**2*VCPP2) ) 

C ***** INVERT [DWDG] ***** 

INV(2 , 1) - 1/ (DWDG( 1 , 2 ) - DWDG(1 , 1)*DWDG(2 , 2)/DWDG(2 , 1) ) 

INV(1 , 1) - - INV ( 2 , 1 ) *DWDG (2,2) /DWDG (2,1) 

INV(2 , 2) - 1/ (DWDG (2 , 2) - DWDG(1 , 2)*DWDG(2 , 1)/DWDG(1 , 1) ) 

INV (1,2) = - INV (2,2) *DWDG (1,2) /DWDG (1,1) 

C ***** COMPUTE NEW RH02 , T2 ***** 

DO 200 I = 1,2 
OLDG(I) - G ( I ) 

DO 210 J = 1,2 

G(I) = G( I) - INV ( I , J ) *W( J ) 

210 CONTINUE 

200 CONTINUE 

RH02 = G(l) 

T2 = G(2) 

ERR = ( (G(l) - OLDG(l) )**2 + (G (2 ) - OLDG(2) )**2)**. 3 
IF (ERR.GT. IE-5) GOTO 150 

C ***** COMPUTE REMAINING THERMODYNAMIC PROPERTIES ***** 

Z2 = 1 + RH02*VB2 + RH02**2*VC2 

P2 - Z2*RH02* (R/MHE) *T2 

U2 - ((PI + RH01*U1**2 - P2) /RH02)** . 5 

CV2 = (R/MHE)/(K-1) - (R/MHE)*T2*(RH02*(2*VBP2 + T2*VBPP2) + 

C RH02**2/2*(2*VCP2 + T2*VCPP2)) 

DPDT = RH02*(R/MHE)*(1 + RH02*(VB2 + T2*VBP2) + RH02**2*(VC2 + 
C T2*VCP2) ) 

DPDRHO = ( R/MHE ) *T2* ( 1 + 2*RH02*VB2 + 3*RH02**2*VC2) 

CP2 - CV2 + (T2/RH02**2)*DPDT**2/DPDRH0 
K2 = CP2/CV2 

A2 = ( CP2* ( DPDRHO )/CV2)**. 5 
MACH2 = U2/A2 
H2 = HT2 - U2**2/2 
CALL VISCOS(T2 ,MU2) 

RE2 = RH02*U2/MU2 
Q2 = RH02*U2**2/2 

S2 = (R/MHE)*(l/(K-l)*LOG(T2) - L0G(RH02) - RH02*(VB2 + T2* 

C VBP2) - RH02**2/2*(VC2 + T2*VCP2)) + SREF 

C = 

c “ - POST -SHOCK STAGNATION CONDITIONS - 
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c - ASSUME POST -SHOCK STAGNATION POINT GAS 

C = BEHAVIOR IS REAL, ENTROPY AND TOTAL 

q = ENTHALPY ARE CONSERVED FROM STATIC VALUES. 

C 

C ================================== 

c **************************************************** 

C * ITERATE FOR RHOT2 , TT2 TO SATISFY CONSERVATION * 

C * OF ENTROPY AND TOTAL ENTHALPY USING TWO-VARIABLE * 

C * (RHOT2 , TT2) SHOOTING TECHNIQUE. INITIAL * 

C * ESTIMATE FROM PERFECT GAS RELATIONS. * 

„ **************************************************** 


TT2 = TO 

RHOT2 = RH02*(1 + (K-1)/2*MACH2**2)**(1/(K-1) )*2 


G(l) = RHOT2 
G(2) = TT2 


350 CONTINUE „ 

CALL VIRIAL (TT2 ,VBT2 , VBPT2 , VBPPT2 , VB3PT2 , VCT2 , VCPT2 , VCPPT2 , 

C VC3PT2) 

W(l) = HT2 - (R/MHE)*TT2*(K/(K-1) + RHOT2*(VBT2 - TT2*VBPT2) + 

C RHOT2**2/2*(2*VCT2 - TT2*VCPT2)) 

W(2) = S2 - (R/MHE)*(l/(K-l)*LOG(TT2) - LOG(RHOT2) - RHOT2* 

C (VBT2 + TT2*VBPT2) - RHOT2**2/2*(VCT2 + TT2*VCPT2) ) - SREF 


RHOT2* ( 2*VCT2 - 


***** COMPUTE JACOBIAN OF MATRIX [W] ***** 

DWDG (1,1) = - (R/MHE)*TT2*( (VBT2 - TT2*VBPT2) + 

TT2*VCPT2) ) 

DWDG (1,2) = - (R/MHE)*(K/ (K-l) + RHOT2*(VBT2 - VBPT2*TT2 - TT2 
**2*VBPPT2) + RHOT2**2/2* ( 2*VCT2 -TT2**2*VCPPT2 ) ) 
DWDG (2,1) = - (R/MHE)*( - 1/RHOT2 - (VBT2 + TT2*VBPT2) - RHOT2* 
(VCT2 + TT2*VCPT2) ) 

DWDG (2,2) = - (R/MHE)*(1/ ( (K-1)*TT2) - RHOT2*(2*VBPT2 

+ TT2*VBPPT2) - RHOT2**2/2*(2*VCPT2 + TT2*VCPPT2)) 


C *** INVERT [DWDG] *** 

INV(2 , 1) = 1/ ( DWDG (1,2) - DWDG(1 , 1)*DWDG(2 , 2)/DWDG(2 , 1) ) 
INV(1 , 1) = - INV ( 2 , 1 ) *DWDG (2,2) /DWDG (2,1) 

INV(2 , 2) = 1/ ( DWDG (2,2) - DWDG(1 , 2)*DWDG(2, 1)/DWDG(1 ,1) ) 
INV(1 ', 2) = - INV (2 , 2)*DWDG ( 1 , 2) /DWDG (1,1) 


C 


310 

300 


***** COMPUTE NEW RHOT2 , TT2 ***** 

DO 300 I = 1,2 
OLDG(I) = G(I) 

DO 310 J = 1,2 

G(I) = G ( I ) - INV(I , J)*W(J) 

CONTINUE 
CONTINUE 
RHOT2 = G(l) 

TT2 — G ( 2 ) 

ERR = ( (G(l) - OLDG(l) )**2 + (G(2) - OLDG(2) )**2)** . 5 
IF (ERR.GT. IE- 5) GOTO 350 
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ZT2 - 1 + RHOT2*VBT2 + RHOT2**2*VCT2 
PT2 = ZT2*RHOT2* (R/MHE) *TT2 


C ***** NEW ESTIMATE FOR MACH1 ***** 

FP(L) - PITOT - PT2 
IF ( FLAG . EQ . 1 .AND. L.LT.2) THEN 
L = L + 1 

MACH1 - MACH1 - DELMACH 
GOTO 510 
END IF 

IF ( FLAG . EQ . 1 ) THEN 
OLDMACH = MACH1 
DELFP - (FP(1)-FP(2)) /DELMACH 
MACH1 - MACH1 - FP(2)/DELFP 
ERR = ABS( (MACH1- OLDMACH) /OLDMACH) 

IF (ERR.GT. IE-5) GOTO 500 
END IF 

C ***** COMPUTE REMAINING THERMODYNAMIC PROPERTIES ***** 

CVT2 = (R/MHE) / (K- 1) - (R/MHE)*TT2*(RHOT2*(2*VBPT2 + TT2*VBPPT2) 
C + RHOT2**2/2*(2*VCPT2 + TT2*VCPPT2)) 

DPDT - RHOT2*(R/MHE)*(l + RHOT2*(VBT2+TT2*VBPT2)+RHOT2**2*(VCT2 
C + TT2*VCPT2) ) 

DPDRHO - (R/MHE) *TT2*(1 + 2*RHOT2*VBT2 + 3*RHOT2**2*VCT2 ) 

CPT2 - CVT2 + (TT2/RHOT2**2)*DPDT**2/DPDRHO 
KT2 = CPT2/CVT2 

AT 2 = ( CPT2* ( DPDRHO )/CVT2)**. 5 
CALL VISCOS (TT2 , MUT2 ) 


C ========= = = ============== ================= 

c 

C = - SONIC POINT PROPERTIES - 

C 

C = ASSUME SONIC POINT GAS BEHAVIOR IS REAL, 

C = ENTROPY AND TOTAL ENTHALPY ARE CONSERVED 

C = FROM THEIR SUPPLY RESERVOIR VALUES. 

C 

C ================================== — ======= ======== 

MACHS = 1 
HTS = HO 
SS = SO 

C ***** INITIAL GUESS AT RHOS , TS FROM PERFECT GAS RELATIONS ***** 

RHOS = RHOO/(1 + (K- 1)/2*MACHS**2)**(1/(K- 1) ) 

TS - TO/ ( 1 + (K- 1 )/2*MACHS**2 ) 

CALL REALEXP (PS , TS , RHOS , ZS , MACHS , US , AS ,HS , HTS , SS , RES , MUS , QS ) 
GOTO 800 

750 CONTINUE 

C =============================================== ========== 

C 

C = - PROPERTIES AT ARBITRARY POINT - 
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c 

c - ASSUME ARBITRARY POINT BEHAVIOR IS REAL, 

C - ENTROPY AND TOTAL ENTHALPY ARE CONSERVED = 

C - FROM THEIR SUPPLY RESERVOIR VALUES. INPUT 

C - MACH NUMBER TO SPECIFY POINT 

C 

C _======================================================- 

WRITE (*,755) 

755 FORMAT (IX, 'ENTER MACH NUMBER AT ARBITRARY POINT IN NOZZLE ') 

READ (*,*) MACHX 
HTX = HO 
SX - SO 

C ***** INITIAL GUESS AT RHOX, TX FROM PERFECT GAS RELATIONS ***** 

RHOX = RHOO/(1 + (K-1)/2*MACHX**2)**(1/ (K-l) ) 

TX = T0/(1 + (K- 1)/2*MACHX**2) 

CALL REALEXP (PX , TX , RHOX , ZX , MACHX, UX, AX, HX, HTX, SX, REX , MUX, QX) 


800 

C 

C 

C 

C 

C 

C 

c 


CONTINUE 

IF ( FLAG . NE . 3 ) GOTO 900 


- PERFECT GAS CORRECTION FACTORS - 

COMPUTE CORRECTION FACTORS TO BE APPLIED TO 
PERFECT GAS CALCULATIONS. 


C *** SUPPLY *** 

TOP = HO/CP 

RHOOP = EXP(CV*LOG(TOP)/(R/MHE) - (SO - SREF)/ (R/MHE) ) 
POP = RHOOP* ( R/MHE )*T0P 
AOP - (K*(R/MHE)*T0P)**.5 

C *** SONIC POINT *** 

TSP = HS/CP 

ASP - (K*(R/MHE)*TSP)**.5 

MACHS P - ( (TOP/TSP - 1)*2/(K- 1) )** . 5 

USP = MACHSP*ASP 

PSP = P0P/(1 + (K-1)/2*MACHSP**2)**(K/ (K-l) ) 

RHOSP = PSP/( (R/MHE) *TSP) 

C *** POST- SHOCK STAGNATION *** 

TT2P = HT2/CP 

RHOT2P = EXP (CV*LOG(TT2P)/ (R/MHE) - (S2 - SREF) / (R/MHE ) ) 
PT2P = RHOT2P* (R/MHE) *TT2P 
AT2P = (K*(R/MHE)*TT2P)**.5 

C *** FREESTREAM *** 

MACH1P = MACHl 

400 DELMACH = MACH1P/10000 

CALL ITER(PT2P , POP ,MACH1P , FI) 

CALL ITER( PT2P , POP , MACH1P+DELMACH , F2 ) 


53 


F = ( F2 - FI ) /DELMACH 

OLDMACH = MACH IP 

MACH IP = MACH IP - Fl/F 

ERROR = ABS((MACH1P- OLDMACH) /OLDMACH) 

IF (ERROR. GT. IE- 5) GOTO 400 

PIP - P0P/(1 + (K- 1)/2*MACH1P**2)**(K/(K- 1) ) 

TIP - TOP/ ( 1 + (K-1)/2*MACH1P**2) 

RH01P - P1P/((R/MHE)*T1P) 

A1P - (K*(R/MHE)*T1P)**.5 
U1P - MACH1P*A1P 
HIP = CP*T1P 
CALL VISCOS (TIP ,MU1P) 

REIP - RH01P*U1P/MU1P 
Q1P - RH01P*UlP**2/2 

C *** POST-SHOCK STATIC *** 

MACH 2 P - ( ( (K- 1)*MACH1P**2 + 2)/(2*K*MACHlP**2 - (K-l)))**.5 
T2P - TT2P/ (1 + (K-1)/2*MACH2P**2) 

P2P = PT2P/ (1 + (K- 1)/2*MACH2P**2)**(K/(K-1) ) 

RH02P = P2P/((R/MHE)*T2P) 

H2P - CP*T2P 

A2P = (K*(R/MHE)*T2P)**. 5 
U2P - MACH2P*A2P 
CALL VISCOS (T2P , MU2P) 

RE2P - RH02P*U2P/MU2P 
Q2P = RH02P*U2P**2/2 

900 CONTINUE 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


**************************************************************** 
-4- TIC CD C C A XT PTinTlT xr K Tm-r-w. Tm 


USERS CAN SUPPLY AN OUTPUT ROUTINE 

* APPLICABLE TO THEIR REQUIREMENTS 

* AT THIS POINT. 


* OUTPUT ROUTINE MUST DO THE FOLLOWING: 

* IF FLAG - 0 OR 1 OUTPUT CONDITIONS AT ALL POINTS IN TUNNEL 

* IF FLAG - 2 OUTPUT CONDITIONS AT ARBITRARY POINT 

* IF FLAG = 3 OUTPUT EQUIVALENT PERFECT GAS CORRECTION FACTORS * 
******************************************************^^^^^ 

*** SAMPLE OUTPUT FOR TESTING PROGRAM *** 

***** FILE OUTPUT ***** 

OPEN (UNIT = 21, FILE = 'RUNDATA.DAT') 

WRITE (21,600) 

600 FORMAT ( IX ,/,' SUPPLY CONDITIONS') 

WRITE (21,605) 

605 FORMAT (2X,'P0 (psi) TO (R) RHOO(slug/ft A 3) ZO') 

WRITE (21,700) P0*1 . 4504E-4 , T0*1 . 8 ,RHOO*6 . 8521E-2/3 . 2808**3 , ZO 

IF ( FLAG . EQ . 2 ) THEN 
WRITE (21,645) MACHX 
645 FORMAT (/, IX, ' CONDITIONS AT M 

WRITE (21,650) 

650 FORMAT (2X,'PX (psi) TX (R) RHOX(slug/f t A 3) ZX' ) 

WRITE (21,700) PX*1 . 4504E-4 ,TX*1 . 8 ,RHOX*6 .8521E-2/3. 2808**3 , ZX 


' , F5 . 2 , ' POINT IN NOZZLE') 
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WRITE (21,655) 

655 FORMAT (IX,' UX (ft/s) MUX (slug/ft-s) REX(l/ft) ') 

WRITE (21,710) UX*3. 2808, MUX*6 . 8521E-2/3 . 2808, REX/3 . 2808 
GOTO 999 
END IF 

WRITE (21,610) 

610 FORMAT (/, IX ,' FREESTREAM CONDITIONS ' ) 

WRITE (21,615) MACH1 
615 FORMAT (IX,' Mlnf = \F5.2) 

WRITE (21,620) 

620 FORMAT (2X,'P1 (psi) T1 (R) RH01( slug/ft "3) Zl') 

WRITE (21,700) Pl*1.4504E-4, Tl*l . 8 , RH01*6 .8521E-2/3. 2808**3 , Zl 
WRITE (21,625) 

625 FORMAT (IX,' U1 (ft/s) MU1 (slug/ft-s) REl(l/ft)') 

WRITE (21,710) Ul*3. 2808, MU1*6 . 8521E-2/3 . 2808, RE1/3 . 2808 

WRITE (21,630) 

630 FORMAT (/, IX, 'SONIC POINT CONDITIONS') 

WRITE (21,635) 

635 FORMAT (2X,'P* (psi) T* (R) RHO*(slug/ft A 3) Z*') 

WRITE (21,700) PS*1 . 4504E-4 ,TS*1 . 8 , RH0S*6 . 85 21E- 2/3 . 2808**3 , ZS 
WRITE (21,640) 

640 FORMAT (IX,' U* (ft/s) MU* (slug/ft-s) RE*(l/ft) ') 

WRITE (21,710) US*3. 2808, MUS*6 . 8521E-2/3 . 2808, RES/3. 2808 


WRITE (21,660) 

660 FORMAT (/, IX, ' POST-SHOCK CONDITIONS ' ) 

WRITE (21,665) 

665 FORMAT (2X,'P2 (psi) T2 (R) RH02( slug/ft *3) Z2 ' ) 

WRITE (21,700) P2*1.4504E-4,T2*1. 8, RH02*6 . 8521E-2/3 . 2808**3, Z2 
WRITE (21,667) 

667 FORMAT (IX,' U2 (ft/s) MU2 (slug/ft-s) RE2(l/ft) ') 

WRITE (21,710) U2*3. 2808, MU2*6 . 8521E-2/3 . 2808, RE2/3. 2808 
WRITE (21,670) 

670 FORMAT (2X,'PT2 (psi) TT2 (R) RHOT2 (slug/ft A 3) ZT2 ' ) 

WRITE (21 , 700) PT2*1 . 4504E-4 , TT2*1 . 8 , RHOT2*6 . 8521E-2/3 . 2808**3 , 

C ZT2 

WRITE (21,675) 

675 FORMAT (IX, ' SHOCK RATIOS ',/, 3X, ' PT2/PT1 P2/P1 ' , 

C ' T2/T1 ' ) 

WRITE (21,680) PT2/P0 , P2/P1 , T2/T1 
680 FORMAT (3X , E10 . 4 , 3X , E10 . 4 , 3X, E10 . 4) 

C *** PERFECT GAS CORRECTION FACTORS *** 

IF ( FLAG . EQ . 3 ) THEN 
WRITE (21,689) 

689 FORMAT (/, IX, ' PERFECT GAS CORRECTION FACTORS') 

WRITE (21,690) 

690 FORMAT (IX,' (P0/Pl)p/(P0/Pl) (T0/Tl)p/(T0/Tl) (RHOO' 

C , '/RHOl)p/(RHO0/RHOl) (A0/Al)p/(A0/Al) ' ) 

WRITE (21,720) P0P*P1/(P1P*P0) , T0P*T1/(T1P*T0) , 

C RHOOP*RH01/(RHOO*RHO!P) , A0P*A1/(A1P*A0) 
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WRITE (21,691) 

691 FORMAT (IX,' (P0/PS)p/(P0/PS) (TO/TS)p/(TO/TS) (RHOO' 

C , '/RHOS)p/(RHOO/RHOS) (AO/AS)p/(AO/AS) ' ) 

WRITE (21,720) P0P*PS/(PSP*P0) , T0P*TS/(TSP*T0) , 

C (RHOOP*RHOS) / (RHOSP*RHOO) , AOP*AS/(ASP*AO) 

WRITE (21,692) 

692 FORMAT (IX,' (PO/P2)p/(PO/P2) (T0/T2)p/(T0/T2) (RHOO' 

C , '/RHO2)p/(RHO0/RHO2) (AO/A2)p/(AO/A2) ' ) 

WRITE (21,720) POP*P2/(P2P*PO) , T0P*T2/(T2P*T0) , 

C RHO0P*RHO2/(RHO2P*RHO0) , A0P*A2/(A2P*A0) 

WRITE (21,693) 

693 FORMAT (IX,' (PO/PT2)p/(PO/PT2) (T0/TT2)p/(T0/TT2) (RHOO' 

C , '/RHOT2)p/(RHOO/RHOT2) ') 

WRITE (21,721) POP*PT2/(PT2P*PO) ,TOP*TT2/(TT2P*TO) , 

C (RHOOP*RHOT2)/(RHOT2P*RHOO) 

END IF 

700 FORMAT ( IX , F10 . 5 , 5X , F8 . 4 , 7X, E9 . 4 , 6X, F7 . 5) 

710 FORMAT ( IX, 3X, F7 . 2 , 9X, E9 . 4 , 11X, E9 .4) 

720 FORMAT ( 8X, F6 . 4 , 11X , F6 . 4 , 15X , F6 . 4 , 15X, F6 . 4) 

721 FORMAT (8X, F6 . 4 , 12X, F6 . 4 , 16X, F6 . 4) 

999 END 


C ############################################################################ 

SUBROUTINE VIRIAL(T , VB , VBP , VBPP , VB3P , VC , VCP , VCPP , VC3P) 

C * COMPUTE 2nd, AND 3rd VIRIAL COEFFICIENTS OF * 

C * HELIUM BASED ON THE CURVE FITS OF EMPIRICAL * 

C * DATA AND THEIR FIRST AND SECOND DERIVATIVES * 

C * WITH RESPECT TO TEMPERATRE . * 

C *************************************************** 

IMPLICIT REAL (A-Z) 

COMMON K,MHE,R, SREF 

C ***** 2nd VIRIAL COEFFICIENT ***** 

IF (T.LE.1300) THEN 
LO = -13.4067 
LI - 165.4459 
L2 = -1357.92 
L3 = 5959.061 
L4 = -12340.8 
ELSE 

LO = 1.178236 
LI = -7.57134 
L2 = 5225.701 
L3 = -188923 
L4 - 2460461 
END IF 

VB - (LO + L1*T**( - 1/4 . ) + L2*T**( - 3/4 . ) + 

C L3*T**( - 5/4 . ) + L4*T**(-7/4. ) )/(1000*MHE) 
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no o o no 


VBP - (-Ll/4*T**(-5/4. ) - L2*3/4*T**(-7/4. ) - 

L3*5/4*T**(-9/4. ) - L4*7/4*T**( - 11/4 . ) )/ 

(1000*MHE) 

VBPP = (Ll*5/16*T**(-9/4.) + L2*21/16*T**( -11/4. ) + 

L3*45/16*T**(-13/4. ) + L4*77/16*T**(-15/4. ))/ 
(1000*MHE) 

VB3P = ( -L1*(45/64)*T**( - 13/4 . ) - L2*231/64*T**( -15/4 . ) - 
L3*585/64*T**( -17/4. ) - L4*11155/64*T**(-19/4 . ) )/ 
(1000*MHE) 

C ***** 3rd VIRIAL COEFFICIENT ***** 

JO - -13.7898 
J1 = 139.7339 
J2 = 8114.259 
J3 = -17456.9 

VC = (JO + J1*T**( - 1/4 . ) + J2*T**(-3/4. ) 

C +J3*T**( - 5/4 . ) )/(1000*MHE)**2 

VCP = ( - J 1 /4*T** ( - 5/4 . ) - J2*3/4*T**( -7/4 . ) 

C - J3*5/4*T**( -9/4 . ) )/(1000*MHE)**2 

VCPP = (+J1*5/16*T**( -9/4 . ) + J2*21/16*T**(-ll/4. ) 

C + J3*45/16*T**(-13/4. ) )/ (1000*MHE)**2 

VC3P - ( - Jl*( 45/64 )*T**( - 13/4 . ) - J2*231/64*T**(-15/4. ) - 
C J3*585/64*T**( - 17/4 . ) )/(1000*MHE)**2 

RETURN 

END 


C 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE REALEXP ( PX , TX , RHOX , ZX , MACHX , UX , AX , HX , HTX , SX , REX , 

C MUX , QX) 

IMPLICIT REAL (A-Z) 

INTEGER I , J , L 

DIMENSION W(2) ,G(2) ,DWDG(4,4) ,OLDW(2) ,INV(2,2) ,OLDG(2) ,FP(2) 
C ,HITER(2) ,SITER(2) 

COMMON K,MHE,R, SREF 


- CONDITIONS AT AN ARBITRARY POINT - 

THE FLOW PROPERTIES OF REAL HELIUM AT AN 
ARBITRARY POINT ARE DETERMINED. THE POINT 
IS IDENTIFIED BY PICKING A MACH NUMBER 
(USUALLY M - 1, THE SONIC POINT, OR M = Minf) AND 
DETERMINING FLOW PROPERTIES FROM CONSERVATION 
OF ENTROPY AND TOTAL ENTHALPY. BECAUSE OF 
THE EXTREME SENSITIVITIES OF SOME OF THE 
DERIVATIVES TO DIFFERENT MACH NUMBERS, THE 
DERIVATIVES REQUIRED IN THE NEWTON ITERATION 
ARE EVALUTED NUMERICALLY INSTEAD OF ANALYTICALLY 
TO PREVENT PROGRAM CRASHING AND TO INSURE 
RAPID CONVERGENCE. 
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G( 1) - RHOX 
G( 2) = TX 

CONTINUE 

CALL VIRI AL( TX , VBX , VBPX , VBPPX , VB3PX, VCX , VCPX , VCPPX , VC3PX) 

CVX = (R/MHE)/ (K- 1) - ( R/MH E ) *TX* ( RHOX* ( 2 *VB PX + TX*VBPPX) + 
RHOX**2/2*(2*VCPX + TX*VCPPX) ) 

DCVXDT = - (R/MHE)* (RHOX* (2*VBPX + 4*TX*VBPPX + TX**2*VB3PX) + 
RHOX**2/2*(2*VCPX + 4*TX*VCPPX + TX**2*VC3PX) ) 

DCVXDR = - (R/MHE )*TX*( (2*VBPX + TX*VBPPX) + RHOX*(2*VCPX + 
TX*VCPPX) ) 

DPDT - RHOX*( R/MHE) *(1 + RHOX*(VBX + TX*VBPX) + RHOX**2*(VCX + 
TX*VCPX) ) 

DPDTDT = RHOX* (R/MHE)* (RHOX* (2*VBPX + TX*VBPPX) + RHOX**2* 

( 2*VCPX + TX*VCPPX) ) 

DPDTDR - (R/MHE)*(1 + 2*RHOX*(VBX + TX*VBPX) + 3*RHOX**2* 

(VCX + TX*VCPX) ) 


DPDRHO - ( R/MHE )*TX*(1 + 2*RHOX*VBX + 3*RHOX**2*VCX) 

DPDRDT - (R/MHE)*( 1 + 2*RHOX*(VBX + TX*VBPX) + 3*RHOX**2* 
(VCX + TX*VCPX) ) 

DPDRDR = ( R/MHE )*TX*(2*VBX + 6*RHOX*VCX) 

AX - (DPDRHO + TX*DPDT**2 / ( CVX*RHOX**2 ) ) ** . 5 
DAXDT = (DPDRDT + TX*DPDT**2/(CVX*RHOX**2)*(l/TX + 2*DPDTDT/ 
DPDT - DCVXDT /CVX) ) /(AX*2) 

DAXDR = (DPDRDR + TX*DPDT**2 / ( CVX*RHOX**2 ) * ( 2*DPDTDR/DPDT - 
DCVXDR/CVX - 2/RHOX) ) /(AX*2) 


W(l) = HTX - ( R/MHE ) *TX* ( K/ ( K - 1 ) + RHOX*(VBX - TX*VBPX) + 
RHOX**2/2* ( 2*VCX - TX*VCPX) ) - (MACHX*AX)**2/2 
W(2) = SX - (R/MHE) *(1/ (K- l)*LOG(TX) - LOG(RHOX) - RHOX* 
(VBX + TX*VBPX) - RHOX**2/2*(VCX + TX*VCPX) ) - SREF 

***** COMPUTE JACOBIAN OF MATRIX [W] ***** 

DWDG (1,1) = - (R/MHE) *TX*( (VBX - TX*VBPX) + RHOX*(2*VCX - 
TX*VCPX) ) - MACHX** 2 *AX* DAXDR 
DWDG (1,2) = - (R/MHE)*(K/(K-1) + RHOX*(VBX - VBPX*TX - TX 
**2*VBPPX) + RHOX**2/2*(2*VCX-TX**2*VCPPX) ) - 
MACHX**2*AX*DAXDT 

DWDG (2,1) = - (R/MHE)*(-l/RHOX - (VBX + TX*VBPX) - RHOX* 
(VCX + TX*VCPX) ) 

DWDG (2,2) = - (R/MHE)*(l/( (K-1)*TX) - RHOX*(2*VBPX 

+ TX*VBPPX) - RHOX**2/2*(2*VCPX + TX*VCPPX) ) 
***** INVERT [DWDG] ***** 

INV(2 , 1) = 1/ (DWDG( 1 , 2 ) - DWDG (1,1) *DWDG (2,2) /DWDG (2,1)) 
INV(1 , 1) = - INV(2 , 1)*DWDG(2 , 2) /DWDG (2 , 1) 

INV(2 , 2) = 1/ (DWDG (2 , 2) - DWDG(1 , 2)*DWDG(2 , 1)/DWDG ( 1 , 1) ) 
INV( 1 , 2 ) = -INV(2 , 2)*DWDG(1 , 2) /DWDG(1 , 1) 

***** COMPUTE NEW RHOX, TX ***** 

DO 120 I = 1,2 



OLDG(I) = G(I ) 

DO 130 J = 1,2 

G ( I ) = G(I) - INV ( I , J ) *W( J ) 

130 CONTINUE 

120 CONTINUE 

RHOX = G ( 1 ) 

TX = G(2) 

ERR - ( (G(l) - 0LDG(1) )**2 + (G(2) - OLDG(2) )**2)** . 5 
IF (ERR.GT. IE-4) GOTO 110 

C ***** EVALUATE PROPERTIES AT RHOX.TX ***** 

CALL VIRIAL (TX,VBX,VBPX,VBPPX,VB3PX,VCX,VCPX,VCPPX,VC3PX) 
DPDT = RHOX*(R/MHE)*(l + RH0X*(VBX + TX*VBPX) + RHOX**2* 

C (VCX + TX*VCPX) ) 

DPDRHO = (R/MHE)*TX*(1 + 2*RHOX*VBX + 3*RHOX**2*VCX) 

CVX = (R/MHE) /(K- 1) - (R/MHE)*TX*(RH0X*(2*VBPX + TX*VBPPX) + 
C RHOX**2/2*(2*VCPX + TX*VCPPX) ) 

CPX = CVX + TX*DPDT**2/(RHOX**2*DPDRHO) 

KX = CPX/CVX 

AX = (DPDRHO + TX*DPDT**2 / ( CVX*RHOX**2 ) ) ** . 5 

UX - AX*MACHX 

HX = HTX - UX**2/2 

ZX = 1 + RHOX*VBX + RH0X**2*VCX 

PX = ZX*RHOX*( R/MHE )*TX 

CALL VISCOS (TX , MUX) 

REX = RHOX*UX /MUX 
QX = RHOX*UX**2/2 

RETURN 

END 

C-- - 

SUBROUTINE VISCOS (T, MU) 

C * COMPUTE VISCOSITY FROM CURVE * 

C * FITS OF EXPERIMENTAL DATA. * 

q *****************************•&** 

REAL MU 

IF (T.LE.1.2) THEN 

MU = (2. 1630-2. 6665E1*T + 1 . 2054E2*T**2 - 1.8741E2* 

C t** 3 + 1. 2682E2*T**4 - 3 . 1823El*T**5)*lE-7 

END IF 

IF (T.GT.1.2 .AND. T.LE.3.6) THEN 

MU = (5.02 - 3 . 2241*T + 2.0308*T**2 - 2.2351E-1 
C *T**3)*lE-7 

END IF 

IF (T.GT.3.6 .AND. T.LE.8) THEN 

MU = (-1.5691 + 3 . 4167*T - 1 . 0317E-l*T**2)*lE-7 
END IF 
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IF (T.GT.8) THEN 

MU - 5 . 023E- 7*T** . 647 
END IF 

RETURN 

END 

SUBROUTINE ITER(PITOT , P01 , Ml , F) 

•*■•*•*•*•*••*■* ic-k'k'k 'k'k'k * ■*** ****** * *** 

* USE NEWTON'S METHOD TO * 

* ITERATE FOR MACH NUMBER. * 
**************************** 

IMPLICIT REAL (A-Z) 

COMMON K,MHE,R,SREF 

PT2P01 - (((K+1)*M1**2)/((K-1)*M1**2 + 2) )**(K/(K- 1) ) 
( (K+l)/ (2*K*M1**2 - (K-1)))**(1/(K-1)) 

F = PITOT/POl - PT2P01 

RETURN 

END 
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